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I.  INTRODUCTION 

A  thorough  study  of  prooassss  such  as  astal  forming,  Impact,  penetra- 
tlon  and  adiabatic  shear  banding  requires  Integration,  with  respeot  to 
time,  of  a  coupled  ayatea  of  nonlinear  partial  differential  equations.  For 
the  aodel  representing  one  of  these  pheaoaena  to  be  soaewhat  realistic ,  it 
should  Incorporate  such  effeots  as  strain  hardening,  strain-rate  hardening 
and  theraal  softening.  These  effeots  are  exhibited  by  aost  aetals  under¬ 
going  large  deformations  at  high  strain  rates.  For  homogeneous  and  simple 
shearing  deformations  of  suoh  visooplastlo  materials,  the  adiabatic  shear 
stress-shear  strain  curve  Is  generally  ooncave  towards  the  origin  and  has  a 
peak  in  It.  At  this  peak,  the  effeot  of  thermal  softening  equals  the  com¬ 
bined  erfect  of  strain  and  strain-rate  hardening.  Under  further  loading, 
the  thermal  softening  overtakes  the  strain  and  strain-rate  hardening,  and 
consequently  the  shear  stress  required  to  maintain  simple  shearing  deforma¬ 
tions  of  the  body  decreases  with  an  Increase  in  shear  strain. 

In  essentially  all  practical  problems  enuaerated  above,  one  needs  to 
integrate  the  governing  equations  well  beyond  the  peak  in  the  stress-strain 
ourve.  Whereas  It  Is  a  trivial  matter  to  carry  out  this  Integration  when 
the  deformations  are  homogeneous,  it  Is  a  rather  tlae  oonaualng  endeavor  to 
do  so  for  non-hoaogeneous  deformations  even  when  the  deformations  are  one- 
dlaenslonal.  Herein  we  dlsouss  our  experience  with  two  methods,  the 
forward-difference  soheme  and  the  Crank-Nloolson  method.  In  eaoh  case,  the 
governing  partial  differential  equations  were  first  reduced  to  a  set  of 
ordinary  differential  equations  by  using  the  Salerkin  finite  eleaent 
method.  Also  in  the  oase  of  the  Crank-Nloolson  method  the  number  or 
unknowns  at  eaoh  point  was  Increased  from  five  to  eight  so  that  only  Tlrst 
order  spatial  derivatives  of  the  unknowns  appeared  In  the  equations.  We 
should  point  out  that  the  governing  equations  are  stifr  and  no  artificial 
viscosity  was  introduced  in  either  oase.  Of  course,  the  Crank-Nicolson 
method  has  artificial  vlsaoslty  Inherently  built  into  It. 

Our  numerical  experiments  reveal  that  the  Crank-Nleolson-Oalerkln- 
Finlte-Kleaent  (CNOFS)  method  allows  the  use  of  tlae  steps  at  least  two 
orders  of  magnitude  larger  than  those  permitted  by  the  Forward-Dlfference- 
Oalerkln-Flnlte-Element  (FOOFN)  soheme  and  still  gives  an  acceptable  stable 
solution.  It  Is  oonoelvable  that  the  efficiency  of  the  forward-difference 
soheae  used  herein  would  improve  If  auxiliary  variables  were  Introduced,  as 
was  done  for  the  Crank-Nloolson  method,  so  that  only  first  order  spatial 
derivatives  appeared  in  the  governing  equations. 

We  refer  the  reader  to  excellent  books1'4  and  references  given  therein 
for  a  discussion  of  various  nuaerloal  integration  techniques.  We  note  that 
Chandra  and  Mukherjee5  reoently  used  the  forward-difference  method  to  inte¬ 
grate  a  stiff  sat  of  partial  differential  equations  soaewhat  akin  to  ours. 
They  used  an  luler  type  soheae  with  autoaatio  time-step  control.  However, 
aeleotlag  parameters  that  oontrol  the  t lme-lnoreaent  automatically  Is  a 
hard  task. 

We  add  that  in  an  earlier  paper6  the  emphasis  was  on  reporting  the 
complete  sat  of  solutions,  obtained  by  using  the  CWQFN  method,  to  equations 
studied  herein.  In  this  paper,  we  provide  detail*  of  the  two  numerical 
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techniques  and  oompare  results ,  for  ona  variable  only,  obtained  by  ualnf 
the  two  methods. 

k .  Formulation  of  the  Staple  Shearing  Problem 

We  study  the  staple  shearing  deformations  of  a  dipolar  vlsoo-plastlo 
aaterlal  and  assume  that  all  of  the  variables  have  been  non-d laenslonal- 
ised.  thus  the  body  oooupies  the  infinite  slab  bounded  by  the  planes 
y  ■  41.  Referring  the  reader  to  Reference  6  for  details,  we  note  that  the 


governing  equations  are 

^  1  “  ly  (  O) 

$  «  ke,yy  ♦  A(s2  ♦  a2)  ,  (2) 

a  «  i»(v,y  -  As)  ,  (3) 

O  •  Mv,yy  ”  J  0)  .  (*»> 


with  boundary  oondltlona 
v(t 1,t)  ■  t  1  , 

6 ,y(± 1,t)  .  0  , 
o (♦ 1 ,t)  •  0  , 


(5) 

(6) 


(7) 

(8) 
(9) 


and  a  suitable  set  of  Initial  oondltlona.  Equations  (1)  and  (?)  express, 
respectively,  the  balanoe  of  linear  momentum  and  internal  energy.  Herein  v 
is  the  veloolty  of  a  aaterlal  particle,  o  its  aass  density,  u  its  shear 
modulus,  (  a  oharaoterlstio  aaterlal  length,  k  its  thermal  conductivity,  e 
its  temperature  ohange  from  that  in  the  rererenoe  configuration,  and  s  and 
a  nay  be  Interpreted  as  the  shear  stress  and  the  dipolar  shear  stress,  k 
superimposed  dot  lndioates  material  time  differentiation  and  a  comma 
followed  by  y  signifies  partial  differentiation  with  respect  to  y.  The 
constitutive  relations  (3)  -  (6)  give  one  possible  model  of  vlsooplastlo 
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materials.  Equation  (6)  Implies  that  the  plastic  parts,  As  and  As/i,  of 
the  strain  rate  and  the  dipolar  strain-rate  vanish  when 

(s2  ♦  a2)*4  <  (1  ♦  *-)n(1  -  a0) . 

*0 

Beoause  of  the  non-dioensional  variables  being  used,  the  initial  yield 
stress  equals  one  in  an  isothermal  and  quasistatlo  referenoe  test.  The 
material  parameters  *  and  n  desorlbe  the  strain  hardening  of  the  material, 
a  the  thermal  softening,  and  b  and  o  the  strain-rate  sensitivity  of  the 
material. 


He  presume  that  the  initial  values  of  9,  a  and  *  are  symmetric  and  of  v 
and  o  antisymmetric  in  y  and  seek  solutions  of  equations  (1)  through  (6) 
with  the  same  symmetry.  Thus  the  problem  is  to  be  studied  over  the  spatial 
domain  [o,l]  and  the  boundary  conditions  beoome 


v(1,t)  >  1, 


v(0,t)  a  0, 


e.yO,t)  «  o, 


e,y(o,t)  >  o, 


o(1,t)  •  0,  o(0, t)  >  o. 

For  the  initial  conditions  we  take 


(10) 

(ID 

(12) 


v(y,o)  >  y,  o(y,o)  ■  0,  *(y,0)  »  i>, 

6(y,0)  a  eQ  ♦  0(y) , 


*(y.o) 


a.  1  (1  ♦*-)  (1  ♦  bAsft)"(1 
0  0 


ae(y,0) ) . 


(13) 


The  values  of  0Q,  sQ  and  *  are  suoh  that,  during  homogeneous  deformations 

of  the  blook,  the  shear  stress  aQ  and  the  strain  corresponding  to  lie  on 

the  shear  stress-shear  strain  curve  for  the  material.  A  in  Eq.  (13)C  is 

given  by  lq.  (4)  with  6  ■  eQ,  a  a  sQ,  +  <*,  o  a  0.  The  function  0 

describee  the  aberration  in  the  initial  temperature  distribution  and  will 
result  in  non- homogeneous  deformations  of  the  body. 


B .  Burner leal  Integration  of  governing  Equations 

1.  Crank-Nloolson-Oalerkln-Pinite-Element  Method. 


With  the  auxiliary  variables 
U  *  V  ,  y  ,  g  a  e,y  ,  P  «  O.y  , 
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we  a an  rewrite  equations  (1)  -  (4)  as 


~r  (s  -  tp),  , 

(15) 

P  ’y 

kg,y  ♦  A(s2  «■  o2)  , 

(16) 

p(u  -  As)  , 

(17) 

lu(u,y  -  j  o)  . 

(18) 

Thus  only  first  order  spatial  derivatives  of  the  unknowns  v,0,s,a,u,g  and  p 

appear  In  he  governing  equations.  Let  H1  denote  the  space  of  funotions 
defined  on  [o,l]  the  square  of  whose  first  order  derivative  is  lntegrable 
over  [0, 1] .  We  approximate  the  unknown  functions  v,0,s  etc,  by  a  linear 
combination  of  the  finite  element  basis  functions  {  $. (y),  1x2,... N}  in  an 

i  1 

N-dlmensional  subspaoe  of  H  .  for  example, 


v(y,t)  >  v^tUjfy).  (19) 

Throughout  this  artiole,  a  repeated  index  implies  summation  over  the  range 
of  the  index.  Using  Salerkin's7  method  we  thus  reduce  equations  (14) 
through  (18)  to  the  following  set  of  equations. 


MijUi 


-QlJVi  • 


(20) 


Mljpl  "  -Qlj°i  * 

HljVl  B  *  l  ^ijpi  ’ 

M1J ®1  *  “VQljsi  *  AiPlj  ’ 

M1J*1  *  u  MlJui  ’  M  Ai*k  Hijk  • 
M1J°1  “  "ut  ^ijul  "  u  Ai°k  Rijk  ’ 


(21) 

(22) 

(23) 

(24) 

(25) 

(26) 


where 
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■a 


(27) 


Mu  f0  W  -  MJi  * 

QU  ~f0  dy  • 

\i  =  Qij  "  (*iV|0  ’ 

Rijk  -  JQ  Vj*k  dy  =  Rikj  =  Rkij  • 

pij  ;//iVs2  +  °2)dy  =  PJi  • 


(28) 

(29) 

(30) 

(3D 


W®  note  that  because  of  the  nonlinear  dependence  of  and  A  upon  s,  o,  ij> 

and  0,  the  coupled  set  of  ordinary  differential  equations  (20)-(26)  is  not 
that  easy  to  integrate.  The  matrices  M^,  Q^,  Q^,  R1Jk  and  P  have  been 

evaluated  by  using  the  linear  basis  functions.  Also  v^(t)  denotes  the 

velocity  of  node  i  at  time  t. 

In  the  Crank-Nicolson  method,  equations  (20)-(26),  assumed  to  hold  at 
time  (t  +  At/2),  are  used  to  preduct  the  values  of  v,0 ,s,o ,g,p,u  and  <f>  at 
time  (t  ♦  At)  from  a  knowledge  of  their  values  at  time  t.  This  is  accom¬ 
plished  by  approximating  ^(t  ♦  At/2)  by  (6i(t  +  At)  -  e^tli/At, 

0^(t  ♦  At/2)  by  (0^(t  ♦  At)  ♦  0^(t))/2,  etc .  and  by  first  evaluating  the 

nonlinear  terms  on  the  right  hand  side  of  (2 0)-(26)  at  time  t.  The  result¬ 
ing  system  of  linear  algebraic  equations  is  solved  for  v^t  At)  etc.,  the 

right-hand  side  in  equations  (20)-(26)  is  now  evaluated  at  time  (t  +  At/2) 
and  the  system  of  equations  solved  again  for  v^t  ♦  At)  etc.  This  itera¬ 
tive  process  is  continued  till,  at  each  nodal  point, 

1^1  +  +  if^l  +  1^1  +  +  l^gl  +  I Ap|  +  | Au |  5  e  (32) 

where  subscript  i  has  been  dropped  from  Vj^  etc. ,  Av  denotes  the  difference 

between  the  newly  found  value  of  v  and  that  used  to  compute  the  right-hand 
side  in  (20)-(26),  and  e  is  a  preassigned  small  number.  The  initial  con¬ 
ditions  (13)  were  used  to  find  v^O)  etc. 
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2.  Forward-Difference-Galerkin-Finite-Element  Method. 


In  this  method  the  field  equations  (1)  and  (2)  were  first  oast  Into 
a  weak  form.  Let  $  and  tf*  be  two  smooth  functions  defined  on  [o,l]  such 
that  <j>(0)  =  <j> ( 1 ) »0 .  With  equations  (1)  and  (2)  multiplied  through  by  4»  and  £ 
respectively  and  with  use  of  the  boundary  conditions  (10)— (12),  integration 
by  parts  over  the  Interval  [0,l]  gives 


(33) 


(34) 


Let  the  Interval  fo,ll  be  divided  into  (N-1)  subintervals,  not  necessarily 

0  1 

of  equal  length.  Thus  N  is  the  number  of  nodes  in  the  mesh.  Let  4>^  ,  4>j 

(1=1,2, ...N)  be  the  Hermite  basis  functions7,  and  4>^  (1=1, 2, ...H)  the 

finite  element  basis  functions  introduced  previously  (e.g.  see  Eqn.  (19)). 
We  impose  the  following  approximations  on  v  and  6. 

v(y,t)  a  v1(t)  4>?(y)  +  Y1(t)  4>j(y)  ,  (35) 

e(y.t)  «  0t(t)  q>i(y) .  (36) 

Here  Y^(t)  is  the  value  of  v,  at  the  node  i  at  time  t.  Hermite  basis 

functions  <>?,  can  be  constructed  by  matching  together  element  shape 

functions  $2  an<1  s*-milarly  4>^( y)  can  be  obtained  by  matching  4>^ 

and  In  fche  Galerkln  approximation,  the  same  set  of  basis  functions  are 

used  to  approximate  the  test  functions  41  and  £  as  are  used  for  v  and  6. 
Recaling  that  equations  (33)  and  (34)  must  hold  for  arbitrary  4>  and  £,  we 
arrive  at  the  following  set  of  ordinary  differential  equations. 

M  w  «  -  F  ,  (37) 

H  6  3  -  ♦  W.  (38) 
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Here 


w  =  {▼1»y1»*?»y 


I'M’  2’r2’  N,tN 


e  =  {e  »e  - 


W  . 


-  s  {F1 »F2»  F2M}  s  1  {fJJ*  f(J+1)J*  f(J+3)J}  ’ 

J=1 


(J+1)J 


'  ( J+2) J 


( J+3)J 


A0  0 
3d)  +  0  <b 


My  Myy 

i  i 

+  0  <j> 

1, y  Myy 

o  o 

♦  a  $ 

2, y  2,yy 


1  1 

a  <p  +  a  <p 

2,y  2,yy 


N-1  r 
M  =  E  / 

Jrl  J 


0  0 

1  0 

0  0 

1  0 

♦ 

♦  ♦ 

* 

1  1 

1  1 

2  1 

2  1 

0  1 

1  1 

0  1 

1  1 

'  4> 

4>  t 

♦  <f> 

$  d> 

1  1 

1  1 

2  1 

2  1 

0  0 

1  0 

0  0 

1  0 

,  t 

♦  4> 

♦  $ 

«>  $ 

1  2 

1  2 

2  2 

2  2 

0  1 

1  1 

0  1 

1  1 

4> 

♦  ♦ 

♦  * 

♦  $ 

1  2 

1  2 

2  2 

2  2 

with  similar  definitions  for  H,  T  and  W.  In  the  above  integrations  fij  is 

the  region  oooupied  by  the  Jth  element.  These  Integrals  are  evaluated 
numerically  by  using  the  4-polnt  Gauss  integration  rule.  Explicit  expres¬ 
sions  for  the  matrices  in  Eq.  (38)  are  not  stated  above  since  they  are 
given  in  many  books  on  the  finite  element  method,  e.g.  Becker  et  al.7 


Equations  (37) »  (38),  and  (3)-(5)  are  integrated  with  respect  to  time  t 
by  using  the  simple  forward-difference  method.  The  solution  of  equations 
(37)  and  (38)  gives  nodal  values  of  v,  y,  and  0  at  the  next  step.  From 
these,  values  of  v,y,9,  and  v,yy  at  the  Gauss  points  of  integration  are 

calculated  by  using  the  interpolation  relations  (35)  and  (36).  For  each 
Gauss  point,  Bqs.  (3)-(5)  are  integrated  to  obtain  the  local  values  of  s,o, 

and  <|»  at  the  next  time  step.  Because  the  integration  scheme  is  only  condi¬ 

tionally  stable  in  the  linear  case,  the  time  step  has  to  be  kept  very 
small;  its  value  is  dependent  on  the  grid  size,  the  material  properties, 
and  the  present  deformations  of  the  body. 

C .  Computation  and  Discussion  of  Res_ults 

In  order  to  compute  numerical  results  the  following  values  of  various 
non-dimensional  parameters  that  correspond  to  a  typical  hard  steel  were 
chosen. 

p  =  3.923  x  10“5  ,  k  =  3.973  x  1 0-3  ,  a  =  0.4973  ,  v  =  240.3  , 

n  =  0.09  ,  <i<0  =  0.017  ,  b  =  5  x  1 06 ,  m  =  0.025  . 

For  homogeneous  deformations  of  the  block,  the  peak  in  the  shear  stress- 
shear  strain  curve  occurs  at  a  strain  of  0.093.  The  uniform  temperature 
0q  =  .1033  in  the  block  when  y  =  0.0692  was  perturbed  by  adding  a  smooth 
temperature  bump 


9(y)  =  0.1  (1-y2)9  e‘5y 

and  the  resulting  initial-boundary  value  problem  wa3  solved  by  using  the 
aforementioned  two  methods.  In  each  case  no  attempt  was  made  to  use 
diagonal  matrioea  equivalent,  in  some  sense,  to  those  computed  by  using  the 
basis  functions.  The  domain  [o,l]  was  divided  into  13  subdomains  with 
nodes  at  0,  .05,  .10,  .15,  .20,  .25,  .34375,  .43750,  .53120,  .6250,  .71375, 
.81250,  1.0.  For  the  forward-difference  scheme  various  integrals  appearing 
in  the  expressions  for  F,  M,  H,  T  and  W  were  evaluated  by  using  the  4-point 
Gauss  quadrature  rule. 

When  t  s  0.0  and  0.01,  the  forward-difference  scheme  necessitated  tak- 

_7 

ing  At  s  .5  x  10  in  order  to  obtain  a  stable  solution.  However,  for  the 

_4 

Crank-Nicolson  method,  At  s  .1  x  10  was  found  to  give  a  stable  and 

-5 

acceptable  solution  since  the  results  obtained  with  At  s  .5  x  10  were 

found  to  be  Indistinguishable  from  those  computed  with  the  larger  value  of 
At.  As  is  clear  from  the  two  sets  of  results  shown  in  Figs.  1  and  2,  the 
non-physical  damping  introduced  by  the  Crank-Nicolson  method  results  in  the 
delayed  response  as  compared  to  that  obtained  with  the  forward-difference 
method.  As  is  dlsoussed  in  Reference  6,  the  development  of  a  late  stage 
plateau  is  a  numerical  artifact  and  does  not  represent  a  physical  phenome¬ 
non.  The  plateau  was  also  developed  in  the  solution  computed  by  using  the 
forward-difference  method  even  though  it  is  not  depicted  in  the  figure. 
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Figure  1:  Comparison  of  Solutions  (for  l  =  0.0)  by  the  Two  Different  Integration  Techniques. 


Figure  2:  Comparison  of  Solutions  (for  l  =  0.01)  by  the  Two  Different  Integration  Techniques. 


The  spatial  variation  of  s  at  late  times  also  indicates  some  kind  of  numer¬ 
ical  instability.  Since  the  average  applied  strain  rate  is  unity,  the 
abscissa  also  represents  the  time  measured  from  the  instant  (denoted  by  I 
in  Pig.  1)  the  uniform  temperature  field  is  perturbed.  On  an  IBM  4381  com¬ 
puter,  the  CPU  time  required  to  compute  the  solution  by  the  finite-differ¬ 
ence  method  was  nearly  three  times  that  needed  for  the  other  method  when  e 
in  Eqn.  (32)  was  set  equal  to  .01. 

Figures  3  and  4  compare  the  solutions  for  t  =  0.0  and  i  =  0.01  obtained 
by  the  Crank-Nleolson  method  and  the  forward-difference  method.  In  each 
case,  g,  =  0.01  results  in  a  delayed  response  in  the  sense  that  yp(0,t) 

begins  to  rise  to  its  maximum  value  slower  and  later.  However,  the  two 
integration  techniques  depict  a  similar  qualitative  difference  between  the 
solutions  of  governing  equations  for  i  =  0.0  and  l  =  0.01.  We  have  plotted 
only  versus  t  in  all  of  the  figures  since  the  yp(y,t)  is  maximum  at 

y=0  and  the  rate  at  which  yp(o,t)  builds  up  is  important  in  physical  pro¬ 
blems.  The  evolution  in  time  of  other  variables,  the  spatial  variation  of 
these  variables  at  different  times,  as  well  as  the  effect  of  choosing  dif¬ 
ferent  perturbations  0(y)  have  been  given  in  [6,8,9]. 

Whether  or  not  the  introduction  of  auxiliary  variables  in  the  FDGFE 
metehod  will  permit  the  use  of  a  larger  time  step  remains  to  be  seen. 

Also,  the  use  of  automatic  time-step  control  as  discussed  by  Chandra  and 
Mukherjee5  may  improve  the  efficiency  of  the  FDGFE  method.  Further  work  in 
resolving  some  of  the  issues  raised  herein  and  selecting  an  optimum  value 
of  At  is  currently  under  progress  and  will  be  reported  on  in  future. 
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Figure  3:  Comparison  of  Solutions  for  i  w  0.0  and  l  *  0.01  by  the  CNGFE  Method. 


Comparison  of  Solutions  for  i  ■  00  and  I  ■  0.01  by  the  POGFt  Nathod 
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non  r\j  noon  noon  non  no  noon 


PR  3GRAU  AO  I  AO  ( INPUT,  JUTP'JT,  TAPES  •INPUT,  'AP£6»0UT  PUT  ) 

dime'isidi  rcaRn<3i),Moo£<2,3i),  enass<io),anass(4,62) 

♦  ,FdRCE!62),EF0*CE<4),  TG00TI31) 

DIMENSION  ID(4)»£TAU(5)»YY(2)»ESIG(5),  MRKf5»)0) 

♦  » VELDDT I  31) ,  DDOOT  1 3, 30) »  SIGMA(5,30) 

♦  ,  PGD1T! 5,  30  ),P')OOT(5,  30 1,  EGD  TT*  3,  30), TEMPI  5,  30) 

♦  , PG AMA I 3, 100 )»SI(5,130), TIP (31),P0(5, 100),TAU(3, 30) 


READ  THE  INPUT  DATA 

REAO<5,1000)  CM, CN, BET A, A,  8,  S  I0,RH  3,  C  V 
»EADC  3,1010)  NT, NtNT,  DT,CK,C J# GO 1T0» C MU, CNU 
READ !5,  1)20)  NUMEL » NODES » NPR IMT, Cl »HT 
REA0!5,132M  TIME#  RST APT 
CMtl-C  1U*1000. 


1000  F  ORMAT ( 8F 10. 4 ) 

1010  FORHATIf’,  13,6=10.4) 

1021  CDPMAT(3I5,2F10.5) 

1023  F0RMATIE10.4*e10.4l 

PRINT  3'r  THE  INPUT  0A*A. 

WR I Tt< 6,2000)  CM,CN, AS TA, A , 1, S 1 0, RHO, C V 
WRITE!  6,  2010)  MT,  HINT,  DT,CK, C J»GD0T0»CMU»CMU 
WRITE  (6,  2020)  IUNEL,NOO£S,CL,HT,NPRINT 
WRITE! 6, ’0  25 )  TINE.RSTART 


GENERATE  MODE  NUMBERS  A'lO  NON-DIMENSIONAL  COORDINATE  5  • 
CALL  GRID! VCORI.N  IDE,  NODES,  NllMEl  ) 


PRINT  nut  the  MON-OINCNSIVIAl  MODAL  COORDINATES. 

00  33  J-l, NODES 
35  WRITE! 6, 2030)  J,TC0RD(J) 

00  40  I«l,  N'JME l 

40  WRITE!  6,2040)  I, NOOEI 1, 1 1,  MODE (2, I ) 

2000  F0PNATI5X,  , E  15.5/ JX,  'N  •'  ,  El  5 , 5 25  K,  •  BETA-*  ,  El 5. 5/ 

♦  5X,'A-',E15.5/3X, *B-',£13.3/3X, 'SI0-'»E15.S/ 

♦  3X» 'RHf)-»,ei5.5/3X, *CV«',  E15.5) 

2110  FORMAT (5T, 'NO,  0*  TIME  STEPS  -*,110/ 

♦  5X*  •  NO,  OF  INTE1.  P1INTS  USED  IM  N'JNERICAl  INTEGRATION  •*,13/ 

♦  5X»  ' TINE  IMCP.ENEMT  •»,E15.5/ 

♦  5X, 'THERMAL  CONDUCTIVITY  .«,E13.5/3X 

♦  , 'FACTOR  TO  CONVERT  rr-)n  JOULES  TO  KG-H  -*,E15.3/ 

♦  3X,  '  •’RESCRIBED  STRAIN  RATE  •  •  ,  El  3 . 3/BX,  *  MU  •  ',£13.3/ 

♦  5X,'NU  -  ',£15.5) 

2120  FORMAT {BX, » NUMBER  OF  ELEMENTS  ««,U0P5X, 

♦  'NUMBER  OF  NODES  •',113 

♦  /5X, ' MATERIAL  LENGTH  ■  ',E15.3/5X, 

♦  'HEIGHT  OF  TIE  SPECIMEN  IH  IN  METERS  )  •  ',  =  15.5 

♦  /5 X,  'PRINT  INTERVAL  - • , T 5 ) 

0  25  FORMAT  (5X»  'TIME  AT  THE  START  0*  rms  R!JN  •',I13.5^3X, 

♦•RESTART  JOB  !=  PSTART  -  3.0;OTh£RWI SE  NOT  A', IX, 

♦  •RESTART  JOB , ' /5 X, 'RST ART  .  «,F13.JI 

COMPUTE  NON-DIMENSIONAL  NUMBERS 
CKAPA.7  •  BF’A*<S  IO»*FNI 


C001C1 
COOIIO 
GOO 120 
000130 
000140 
000150 
000160 
000170 

ooo  ieo 

000140 
COO  200 
000210 
C00220 
C00239 
000231 
000240 
000251 
0003CO 
000330 
000340 
000350 
0C03S1 
ooo  3  eo 
000340 
C0O4C0 
000410 
000420 
000430 
000411 
000440 
000450 
000460 
000470 
000400 
000440 
000500 
COO  5 10 
OOO  3  20 
0003)0 
000340 
000510 
000560 
000570 
000560 
000540 
000600 
000610 
000620 
0006)0 
000640 
000650 
000660 
000670 
000600 
000640 
000700 
000710 
000711 
000712 
000713 
000720 
OOO  t  JO 
000740 
0007  50 


20 


« 


■  »  -  .KM'Mvftr  I  i-«.  »*\.J  I  0007 to 

A  -  A*TR  000770 

PACT!  •  4.0*CK/IHT*HT*RH9*CV*GnST0>  000700 

FAC l  •  2«*CL/'IT  000700 

OHO  •  RHn*UHT*GMTO/2.t**21/CKAFAO  000900 

RH9I-1.0/RH0  000009 

CRU  •  CNU/CKAFAO  000910 

CNU  -  CNU*2./(CKAFAO*HT*CL)  000020 

VR1T*I6,2099)  CK AFAO, F ACT1»RH0,CNU»CHU# A#  TO  000090 

2099  FORMAT ( 9X» *KAFA0  •  *»E19.S/9X, ‘THERMAL  FACTO*  •  *,E19.9/5X,  0000*0 

♦  •!NE”tA  FACTTR  •  *,E19.9/9X,  *NU9AR  ■  *,E19. 9/9X,  000090 

♦  'HU-4AR  •  *»E19.9/9X,  •MO'I-DTNENST'JMAL  A  ■  SE19.9  COO 9 60, 

♦  /9X,  'REFERENCE  TE4f EOATlUE  •  NE14.4)  COOOToi 

2090  FORMAT (5X, 19#  F19.9)  COOOOO 

2044  FQRNATI9X, 91191  000090 1 

C  000400 1 

C  HE  AO  THE  FARASET E*S  COtTRlLLING  THE  DIST'MO  ANCE.  C00910 

C  000420 

R(AO( 9*10001  ALF  A.CNNfEFSILON  000990 

REAOI  9*10901  HT  A'J#  HOF#  HS  I, HGF09T»HT4*/N0T  000491 

C  READ  (9# 1091)  TI NE.HTAU#  HO*#  HOT#  HT4F#HJ1#  H6*90T  000992 

C1411  F9RNATI50X/E29.10/10X,  E20# 10/241# E20# 10/ 9 IX# E 20# 10/  000499 

C  *  29I,E20.1G/14X,E29.10/9JX,E20.10)  000994 

C  VRITEI  4#2091 1  HT  Ai)#HO*#HOT#HTH*»  HSl»HOFDOT  000499 

MRIT2I 4,20901  AL  F  A,CNN,«FX  IL  TN  000440 

2090  FORMAT ( 9X#  'VALUE  9  OF  VARIIDLSS  CONTROLLING  THE  0ISTUR4ANCE  »/  000990 

♦  9X» 'ALFA  •  *,Ei9.9/9X,'N  >  E19. 9/9X, 'EFXIL3N  •  *,119.91  000460 

C  C00970 

C  OEHERATE  VALUES  0*  SHAFE  FUNCTIONS  AND  THEIR  DERIVATIVES.  000909 

C  000440 

CALL  SHARES  INIHTJ  001000 

C  001010 

C  CALCUALTE  VALUES  0*  0 ANA-0 OT  AHO  D-OOT  CAUSED  DT  001020 

C  THE  INITIAL  OXSTUROANCE.  001090 

C  001040 

1F(RSTART,NE.0.0)  001041 

•CALL  9 IST,,FO (TC  ]R0,H<I0E*  EOOOT, DOOOT,  ALFA,CrDI.  EFS1L0N  001090 

♦  »HUNCL,NINT.TCN*.TN».H40CS,HTNF  )  001060 

C  001070 

C  CALCULATE  the  INITIAL  TIELO  stress.  001000 

CK  •  lll.O*HSI/SlO)**CN)*ll.O*0*ODOT0*HO*O*rr)**CM  001040 

C  001100 

C  T ADO  •  T AMO/CKAF 40.  001110 

C  001120 

TAUO  •  HVAU  001190 

C  001140 

C  SET  INITIAL  STRESS  •  TIEL9  STRESS.  THE  INITIAL  OISTUROANCE  ALTERS001 190 

c  the  to*al  strain  rate  at  a  node  foint.  oouto 

C  001170 

C  SE*  the  INITIAL  STRAIN  RATI  «  STRAIN  RAT*  DOE  TO  HONOGENEOUS  OOllEO 

C  DEFIRNATION  *  STRAIN  RAtX  C4USE0  or  THE  OISTUROANCE.  001140 

C  NON-MlNEHSTlHALItn  STR4IN  RATE  DUE  TO  N0N1GSNE9US  OEFOR NATION  0012CO 

C  IS  EOMAL  TO  1.9  001214 

C  001220 

IFIRSTAR*.E4.0.0)  00  TO  99  001221 

90  49  NBL  •  l,  UUNEL  001290 

DO  44  t)(T  •  1,  NIHT  001240 

(009T(II|T,HEL)  •  E«04T«INT,NIL)  ♦  1.0  0C1290 

4*  CONTINUE  001240 

C  SET  INITIAL  STRESS  •  STRESS  CAUSED  DT  THE  DISTUR44MC?  ♦  TAUO  D012T4 

C  STRESSES  CAUSE*)  OT  THE  OlSTURtANCE  ARE  TAKEN  EO*IAL  TQ  0  C012EO 

C  0C124O 

NO  44  t  >  1,NMN(L  001  ICO 

40  4A  2  •  l,HIir  001110 
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f 


r»r>nnr»on^>  o  r»  o  o  n  O  n  n  n  non  o  r>  f*>  O  non 


TR  •  CKAPAO/ (RH3*CV*CJ )  C00760 

*  .  a-TR  000770 

FACU  •  4.0*CK/(HT*HT*RH0*CV*G03T0)  000700 

FACl  •  2.*CL/MT  000790 

RHO  •  RHO*  ( ( HT*0D0T0/2  « )**2)  /CKARAO  000800 

RM0I-1.0/RH0  000003 

CHU  •  CNU/CK ARAD  000310 

C1U  -  C1U»2./ICKARAO*HT*CL)  000820 

WRITS(6*2055)  CKAR40*FACn*RH3,C11U,C1U,A*TR  000830 

2033  F 0R1AT I 9X» 'KARA0  •  •* £13. 3/3X, 'THERMAL  FACTO*  •  '»E19.5/5X»  0C08*0 

♦  •  !M£8TTA  FACT1R  •  '*E19.3/9X» 'HU3AR  •  '*E13.3/9X*  C00830 

♦  •  HU— 3AR  ■  '»E19.3/3X»  'lON-DIMENSTTNAL  A  •  ',E13.5  COOOOO 

♦  /9 X*  • REFCREMCE  TEMPERATURE  •  **E19.9>  CC0070' 

2030  FORMAT ( IX/ 13#  F13.5)  COOOOO 

204.1  F3RHAT(5X*3I10)  000090 

000900 

REAO  the  PARAMETERS  CQHtR3LLIMG  THE  DISTURBANCE.  C00910 

COO  9  20 

RE  ADI  9.1000)  ALFA*CMN,ERSILON  000930 

REAOI  9*13301  HTAU»HGR*HSI*HGi»OOT*HTM**HGT  000931 

READ  (9,1031)  TI ME.HTAU* HGR#HGT»HTMF»HSI*H6*00T  000932 

F IRHAT 1 ?OX*E20.10/10X#  E20. 10/24X#  £20.10/ I3X*  E20* 10/  C00933 

♦  23X*E29.10/14X*E23.10/SJX,F20.101  000934 

VRireiWOll)  HT  A'I,HGR»UGT,HTM**  HSI*HGRDOT  C00939 

WRITE!  4*2030)  AL  FA, CNN,  EPS  XL  Tl  C00940 

FORMAT ( -X, 'VAL^IES  0*  VARIISLES  C  31  TROLLING  THE  DISTURBANCE  •/  0009-0 

♦  3X,'ALFA  -  «,E19.9/5X,'N  -  *»E13. 3/3X, ' EPSILON  •  '*E13.3)  C00960 

C00970 

GENERATE  VALUES  0*  SHARE  FUNCTIONS  ANO  THEIR  DERIVATIVES.  000900 

000990 

CALL  SHARES  (MINT)  OOIOCO 

001010 

CALOJALTE  VALUES  OR  GANA-DOT  ANO  D-OOT  CAUSED  IT  001020 

The  initial  DISTURBANCE.  001030 

001040 

IF(A5T*A”.he.9.0I  001041 

♦  CALL  01  ST' ira (YC 3RD, NODE, E6D0T, DODOT,  ALFA* CNN*  EPSILON  001090 

♦  .M')HEL*HIHT,TE'1*,TN**H0DES*HT1R  )  001069 

001070 

CALCULATE  THE  INITIAL  TIELO  STRESS.  OOIOOO 

CF.  •  I  (1.0,HSI/SIOI**CN1*(1.0*9RODOTO*M«®03T)**CN  001090 

ooneo 

r alio  •  tamo/crar\o.  oomo 

001120 

TAUO  •  H’AU  001130 

001140 

SET  INITIAL  STRESS  •  TIELO  STRESS.  THE  INITIAL  OISTUROANCC  ALTERSOOUSO 

the  TO”AL  STRAIN  RATE  AT  A  NODE  ROINT.  001160 

001170 

5ET  tie  INITIAL  STRAIN  RATE  •  STRAIN  RATE  DUE  TO  HOMOGENEOUS  001100 

•>EF1*HA-I3H  ♦  STRAIN  RAT*  CVISE9  RT  THE  DISTURBANCE.  001190 

NDN-DlrE'IStTNALireO  STRAIN  RATE  ONE  TO  HOMOGENEOUS  DIFORNATIO*  001  ICO 

IS  EOMAL  T3  i.S  001210 

C  001220 

IFIRS’AR’.EO.O.O)  GO  to  99  001221 

NO  43  NEL  •  1,  IUNEL  001230 

00  *9  lr|T  ■  L,  HINT  001240 

EG01T( IHT» NED  •  EG03T<»NT*NEL)  ♦  1.0  0C12S0 

49  CONTINUE  001260 

C  SET  INITIAL  STRESS  •  STRESS  CAUSED  IT  THE  DISTURBANCE  ♦  TAUO  001270 

C  STRESSES  CAUSES  R»  THE  DISTURBANCE  ARE  TAKEN  E4*ULL  TO  0  C01200 

C  0C1 290 

on  4.6  I  ■  1  *  nunc  L  001  ICO 

90  46  J  •  l.  Hill”  001310 
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SI3MAU/I)  ■  0.0  00X320 

46  TAUIJ/I)  -  CK*<1.0-A*TEMP!J/I) )  001320 

55  CONTINUE  00X331 

C  FOR  THE  RESTART  JOB#  READ  VELUCITY/TEM^SI/TAU/S IGNA.  001332 

IF<RSTART.ME,0.0>  GO  TO  TO  00X333 

00  60  NEL  •  1/NUMEL  001334 

00  60  IHT  -  1/MINT  ''01335 

REAO( 5/10?0)  TAIJ(INT/NEL)/SIGNA(INT/NEU/PGAMA(INT/NEU  001336 

READ (5/ 1071)  SIC INT/MEL)z°GDQT!IMT/MEL )/TEMP t INT/NEL )/ PDCINT/NEL )  C01337 
READ!5/1C72)  EGD4T! INT/NEL )/ DDDOTI INT/NEL)  001338 

60  CONTINUE  001330 

1070  FORMAT ( 10X/3E20. 12)  001340 

1071  FORMAT ( 2X/ 4E19.12)  001341 

1072  FORMAT (2X/2E10.1 2)  001342 

DO  65  NOD  •  1/  NODES  001343 

65  REA0(5/1075)  VELDOT(NOO) /TG0QT(NQD )/TMR(NOD )  001344 

1075  F0RMAT(15X/3E15.7)  001345 

CALL  DD(DDD0T/NINT/MIJMEL/M0DE/VELD0T/TGDJT/YC0RD/EGD0T/DT/N2ER0/  001346 

♦NTT)  001347 

70  CONTINUE  001348 

CMI-l./CM  001350 

COOT  31  •  1 .0/GDQT0  001355 

C  001370 

C  001390 

1030  FORMAT ( 6F13.4)  001396 

WRITE (6/ 2031)  NT AU/HGP/HSI/HGPDOT/HTMP/HGT  001397 

2031  F3RMAT(5X/»MTA*J  -*/E13.5/5X/ »HG“  - */£15.5/SX/  •  HSI  -'/E15.5/  001398 

♦  DX/'HGPDOT  ■•/E15.5/3X/ •  HTMP  -•/ E15.5/5X/ • H6T  -'/E15.3)  001399 

C  001460 

NOF*2*NOOES  001470 

MB-4  001480 

00  990  NTT-l/NT  001490 

C  001500 

C  INITIALIZE  T'lE  GLOBAL  MASS  MATRIX  AND  THE  GLOBAL  FORCE  VECTOR.  001510 

C  THE  MASS  MATRIX  IS  GENERALIZED  IM  THE  SENSE  THAT  IT  INCLUDES  001520 

C  THE  INERTIA  TERMS  APPROPRIATE  FOR  D-DOUBLE  DOT.  THE  FORCE  VECTDR001530 

C  EQUALS  THE  RESULTANT  OF  FORCES  CAUSED  BY  STRESSES  ANO  001540 

C  DIPOLAR  STRESSES  DUE  TO  THE  DISTURBANCE.  001550 

C  0C156O 

DO  50  1-1/ HOF  001570 

FORCE! I)«0.0  001580 

IF(NTT.NE.l)  63  7Q  jo  001590 

00  49  J>1/NB  001600 

49  AMASS! J/ I)  •  0.0  001610 

90  CONTINUE  001630 

C  001640 

C  „  001650 

TIME  •  TIME  ♦  DT*G01T0I  001660 

C  001670 

HCOUMT  •  NTT/MPRIMT  001680 

NZER4  •  MCD'BIT  *  NPRINT  -  NTT  001640 

C  RIMO  THE  -HMOGSHENUS  SOLUTINN  FOR  THE  PRESENT  VALUE  OP  TIME.  001 7C0 

C  001710 

CALL  MJ.NT6HTHP/HGT/46P  /  !TA'I,  HSI/CHI/  CM/A/B/GOOTO  001720 

♦  /  CNM/OTzSIO/TINE,  H6PODT/NZERO)  001730 

C  001740 

C  001746 

C  ASSEMBLE  THE  GLOBAL  MASS  MATRIX.  001750 

C  001760 

DO  900  N€LM*  1/NUMEL  001770 

30  104  1*1/2  001780 

II*NQOE! I/NCLM)  001790 

TY!I)  -  YCORD! II )  001 800 


100  CONTINUE 


001810 
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00  105  INT  ■  1*  NIMT  001820 

ETAUdNT)  -  TAU(  INT,NELM)  001830 

105  ESIGdNT)  -  SIGN A (INT, NELM)  001840 

001850 

CALCULATE  THE  ELEMENT  MASS  MATRIX  AND  THE  FORCE  VECTOR  001860 

001870 

CALL  ELEMEMTfYY,  ETA'J,  ESIG,  EFQRCE, E HASS, NIMT, F ACL, NTT )  001880 

no  110  11-1,2  001890 

JJ«2*II  001900 

I  •  M1DE(II,NELM)*2  001910 

LD(JJ-1)-I-1  001920 

110  L0(JJ)-I  C01930 

MC-0  0019 AO 

00  13)  11-1,4  001950 

I-LOdI)  001960 

FORCE(I)  -  FORCE ( I )  ♦  EFORCE(II)  001970 


IF(MTT.ME.l)  GO  TO  130  001980; 


DO  120  JJ-l, II  001990 

MC-MC+1  002000 

M-LD(JJ)  002010 

IR-NIN0d,N)  C02020 

IC-IABS(I-M)*1  002030 

120  AMASSdC/IR)  -AMASS  (IC, IR )  ♦ENASS (HC 1  0020AO 

130  CONTINUE  002060 

900  CONTINUE  002070 

00  915  J-1/M0F  002080 

FORCE t J )  -  FORCE (J)*RH0I  C02085 

915  CONTINUE  002G90 

002110 

MODIFY  THE  NASS  MATRIX  FOR  THE  PRESCRIBED  BOUNDARY  CONDITIONS.  C02120 

002130 

AT  THE  BOUNDARY  NODES  (THE  FIRST  ANO  THE  LAST  NODE)  0021AO 

GAMA  ANO  VELOCITY  IS  PRESCRIBED.  002150 

THEREFORE,  ME  SET  GAMA  -  DOUBLE-DOT  -  0.0  AND  ACCELERATION  -  0.0  002160 

AT  THESE  NODES.  002170 

002180 

II  •  2*N0DES-1  002240 

CALL  HDOIFYt 1,0.0, AMASS, FORCE, MOF, MB, MB )  002250 

CALL  NOD IF Yd I, 0 ,0, AM ASS, FORCE, NDF,N9, MB )  002255 

IF(NTT .EQ.l)  002265 

♦CALL  SOLVE ( AMASS, FORCE, NOF, MB, MB, 1 )  002270 

CALL  SOLVE ( AMASS,FORCE,NDF, MB, MB, 2 )  002280 

002290 

002300 

002310 

NOW  THE  ARRAY  FORCE  CONTAINS  TOTAL  MODAL  ACCELERATIONS  AND  NODAL  002320 


STRAIN  DOUBLE-DOTS  CAUSED  BY  THE  DISTURBANCE.  eROH  THESE  VALUES  FIND002330 
'■VALUES  OF  DDOT.  002340 

002350 

TRANSFER  VALUES  from  THE  ARRAY  FORCE  INT  ARRAYS  VELDOT  AND  ODDOT.  002360 

002370 

2110  *0RMAT(5X, ’ACCELERATION  DUE  tq  DISTURBANCE  AND  TOTAL  GAMA-2  DOT*/  002390 


♦  5X, ’NODE  #'»5X, 'ACCELERATION*, 5X, *GAMA-DOUBLE  DOT*)  C02400 

DO  200  1-1,  HIDES  002410 

II  -  2*1  002420 

VELDOT(I)  -  FORCE(II-l)  002430 

TGDOT(I)  •  FORCEdI)  002440 

2100  PDRMAT(5X, I3,2E15.5)  002460 

200  CONTINUE  002470 

002480 

002510 

FIND  PLASTIC  PARTS  OF  EGDOT  AMO  DDDOT.  DENOTE  THESE  002520 

BY  PGOflT,PDDOT.  002530 

002540 


uuuu 


002550 
002 560 
002541 
002570 
C02571 
002572 
002573 
002620 
002630 
002640 
002650 
002660 
002661 
002670 
002680 
002690 
002691 

WRITE (6/ 2401 )  SI  (INT/NEL >/ PGDOT l INT/ MEL */ TEMP ( INT/NEL ) /PD( INT/NEL* 002692 


WRITE (6/2401*  EGD0T( INT/NEL*/ DDD0T( INT/ MEL)  002693 

985  CONTINUE  C02694 

on  980  nod-I/Mooes  002695 

980  WRITE (6/2290*  tnD/VELDOT{NaO)/TGOOT(N3D*/TMP(N3D*  C02696 

2290  FORNATdOX/  I5/3E15.7*  002697 

2400  FORMAT ( 27./ 214, 3E 20. 12)  002698 

2401  BORMAT( 2X/  4E19.121  C02699 

990  CONTINUE  002730 

STOP  002740 

EMD  002750 

SUBROUTINE  HOMOG(HTMP/ HGT/ HGP  /HTAU/HSI/CMI/  CN/A/9/GD0T0  C02760 

♦  /  CNU/DT/SIO/TIWE/  HGPDOT/NZEROJ  002770 

HCKAPA1  •  (1.0  ♦  HSI/SIO)**CN  002780 

HCKAPA  -  HCKAPA1  •  (1.0  -  A  *  HTMP1  002790 

OT-HTAU/HCKAPA  C02800 

IF (QT.LT.l  .0 )G1  TO  500  002810 

R2-0T**CMI  -  1.0  002830 

GAMA  -  R2/(8*G00T0  *HTA'JJ  002840 

GO  T3  600  002850 

500  WRITE(6/ 2000)  QT/TIME  002860 

GAMA  -  0.0  002870 

600  CONTINUE  002880 

HSIDOT  •  GAMA*(HTAU*HTAU/H,;KAPU)  002890 

,  HGPDOT  -  GAMA  *  HTAU  002900 

HTHPOOT  »  HGPOQT  *  HTAU  C02910 

HTAUOT-CMU*(1.0-HGPOOT)  002920 

HSI«HSI*HSIDOT*OT  002930 

HTAU»HTAU+HTAUOT*OT  002940 

HGP*HGP*HGPOOT*0T  C02950 

HGT-HGT*  DT  002960 

HTMP  -  HTMP  ♦  HTMPOOT9DT  002970 

IF (NZERQ.EQ.O) WRITE(6/2005 )  TIME/ HTAU/ H6P/ HGT/ HTMP/ HSI/HGP COT  002980 

2005  FQRMAT(5X/ 'HOMOGENEOUS  SOLN.  AT  "HYSICAL  TIME  (SECS.*  -'/E 20.10/  002990 

♦  5X/*TAU  -'/E20.10/5X/ »GAUA  -PLASTIC  •'/E20.10/  003000 

♦  5X, 'TOTAL  STRAIN  •'/ E20.10/5X/ 'TEMPERATURE  •  '/E20.10/  003010 

♦  5X/»SI  •  '/E20.10  003020 

♦  /?X/ •  PLASTIC  STRAIM-P.ATE  -'/E20.10*  003021 

RETURN  003030 

2000  F0RNAT(3Y/ 'FOR  THE  HOMOG.  SUN.  THE  MATERIAL  IS  DEFORMING*  003040 

♦  /lX/'ELASTICALLY'Z  003050 

♦3X/ 'RATIO  OF  STRESS  TO  HARDENING  FUNCTION  -»/F13.3/  CC3060 

♦3X/  'CURRENT  VALUE  OF  TIME  ■  SF15.5)  003070 

END  003080 

SUBROUTINE  SHAPE5(NINT*  003090 

CQMM0H/SHAPE5/$HAP0(2/ 5)/SHAPl(2/5)  /  DSHAPO( 2/ 5*# DSHAP1 (2/5  1/  003100 

♦  ODSHAPO (2/51/00 SHAP1 ( 2/  5 )  / SNAP  B ( 2/ 5 )  / DSHAP3( 2/ 5  *  003110 


CALL  PLASTIC (SGOOT/OODOT/ PGD7T/ P0D7T/  C.NI/G03T0/ SIGMA/  CN/ 

♦  A/B/CMII/CNU/DT/TaUO/TEMP/SIO/  NTT/ NUHEL/ HINT/ WORK/ TAU 
♦/HSI/ WP/HGPOOT/ "G AM A/ SI/R START/ PD  * 

C 

C  CALCULATE  TOTAL  DOOT  AND  TOTAL  GDOT. 

CALL  OD(DDD3T/NINT/NUHEL/N30E/VELDOT/TGD3T/YCORD/BGDOT/OT/NZERO/ 

♦  NTT* 


SOLVE  THE  THERMAL  PROBLEM 

CALL  THERHITEMP  /W0RK/FACT1/  NUHEL/NIMTz NTT, DT/TM® 

♦  /VC1R0,  NODES/ -lonE* 

IF(NZERO.  NE.  01  GO  TO  990 
DO  985  MEL  -  1/ilUMEL 
DO  983  IHT  -  1#  NINT 

WRITE (6/24001  MEL# INT/TAIH INT/ NELJ /SIGMA ( INT/NEL )/PGAMA( INT/NEL) 
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DIMENSION  St  (4) 

SKI)  -  -0.861X36311504033 
SI(2)  »  -0.33OO81043584856 
SK3)  -  -  SIC2) 

SK4>  -  -  SKI) 

EVALUATE  THE  SHAPE  FUNCTIONS  AMO  THEIR  DERIVATIVES. 

DO  100  I  -  l/'IIMT 
S«SI(I) 

SHAPOd*!)  -  (2.0  -  3»0*S  ♦  S**3)/4.D 
SHAP )( 2/ I )  -  (2.0  ♦  3 .0*S  -  S**3)/4.0 
SHAP1(1*I)  •  (1.0  -  S  -  S*S  ♦  S**3 )/4,0 

SHAP1( 2*  T)  -  (-1.0-  S  ♦  S*S  ♦  S**3)/4«0 

DSHAP0(1/I)  •  (-3.0  *3.0*S*S)/4.0 
DSHAP3(2*I)-  (3.0  -  3.0*S*S)/4.0 
0SHAP1 (1* I )-  (-1 .0-2. 0*3*3. 0*S*S ) /4.0 
DSHAP1(2*I>-  (-1.0  >2.0*5  *3.0*S*S)/4.0 
DDSHAPOd* I)  -  (6.0*S)/4.0 
00$HAPO{2» I)  -  (-6.0*S )/4. 0 
DDSHAP1 (1/ I )  »  (-2.0*6. 0*S)/4.0 
DOSHAP 1( 2*1)  -  (  2.0  ♦  6.0*S)/4.0 

SHAPB(1#I)  -  (1.0  -  S)/2.0 

SHAP9(2*I)  -  (1.0  ♦  S)/2.0 

DSHAPB(1*I)  •  -1./2. 

DSHAPB (2*  I )  -  1./2. 

C 

100  CONTINUE 
RETURN 
END 

SUBROUTINE  ELEMENT (YY, ET  V),ESIG*EF0RCE»EHAS5*nlN7*FACL*NTT ) 
COHNOM /SHAPE  S/S HAP 0(2#  5)*SHAP1(2*5)#0SHAP0(2*5)#  DSHAP1 (2*5)# 

♦  DDSHAP0(2*5)*  DDSHA°1 (2*5)*SHAPB(2*5)*DSHAPB(2»5) 

DINENSIOM  YY(2)#  ETAU ( 5)*  V)E IGHT (5 )*  EFQRCE( 4)#EMASS(10l»ESIG(5) 
WEIGHTd)  -  0.34-»85484513''454 
WEIGH* (2 )  -  0.632145154362546 
WEIGHT (3)  -  WEIGHT(2) 

WEIGHT(4)  •  WEIGHT (1 ) 

MC-0 

C 

DO  100  I  -  1#  4 
EFORCE(I)  -  0.0 
IF(NTT.ME.l)  GO  TO  100 
00  50  J-  1*  I 
MC«MC*1 

30  ENASS( MC ) *0.0 
100  CONTINUE 

DO  200  IMT  -  1*N INT 
WT-WEIGHT(IHT) 

DJ AC  -  (YY(2)  -  YY(1) )*0.5Q 
C  WRITE (6# 2010)  DJ  AC 
2010  FORMAT ( 5X#  *  OJAC  -  *#E15.5) 

IF(NTt.he.I)  60  TO  185 
ST2  •  WT*DJAC 

ENASS(l)  •  EMASS (1)  *  SHAPO( 1# IMT) *SHAP0( 1* IMT )*ST2 
EMASS (2)  •  EMASS (2)  ♦  SHAPO (1, IMT) *SHAP1 ( 1# INT )* ST2 
EMA$S( 3)  •  EMASS (3 )  ♦  SHA»1(1» INT)*SHAP1(1# IMT)*ST2 
EMASS(4)  -  EMASS (4)  ♦  SHA»0(1#INT)*SHAP0(2# INT)*ST2 
i;  MASS (  3 )  -  EMASS  (3)  ♦  SHARK  1*  INT)  *SHAP0(2»  1ST  )*ST2 
EMASS! 6)  •  EMASS (6)  ♦  SHAP0(2# IMTI*SHAPO( 2* INT)*ST2 
ENASSf?)  -  EMASS (T)  ♦  SHAP0(1» INT)*SHAP1(2# IMT)*ST2 
EMASS!  9)  •  EMASS  (8)  ♦  SHAPld*  IMT)  *SHAP1(  2*  IMT  )*ST2 
EMASS! 9)  •  EMASS!*)  ♦  SH4»0(2*INTJ*SHAP1(2*INT)*ST2 


003130 

003160 

003170 

003100 

003ZC0 

003210 

003220 

003230 

0032AO 

003250 

003260 

003270 

003250 

Q03290! 

0033C0 i 

003310 ' 

003320 

003330 

0033AO 

003350 

003360 

003370 

003380 

003390 

003400 

C03410 

003420 

003430 

C03590 

003620 

003630 

003640 

003650 

003660 

003680 

003710 

003720 

003740 

003750 

003800 

003810 

003820 

003830 

003840 

003950 

003860 

003870 

003880 

003890 

0039C0 

003910 

003920 

003930 

003940 

003950 

003960 

003970 

003980 

003990 

0040C0 

004010 

004020 

004030 

004040 


non 


ENASSC10) »  EMASS (10)  +  S4A»1 t Z, INTI *SHAP1(2» I(|T)*ST2  00*050 

185  CQNTIflUE  00*060 

IF(OJAC.LE.l.OE-’O)  GO  TO  3000  C0*070 

WT1«WT*FACL/D4 AC  00*080 

GO  TO  1<>0  00*090 

3000  WRITE! 6* 2020)  OJAC  00*100 

2020  FORMAT (5X/ *DJAC  •  SE15.5J  00*110 

STOP  00*120 

1<»0  CONTINUE  00*130 

EFORCE(l)  -  EFORCE(l)  -  (ETAU( INT) ROSHAPOC 1> I NT) *WT  ♦  00*1*0 

♦  esigcint)*dds!upo<i,int)*vti)  co*i;o 

E FORCE  12)  ■  EFORCE (2)  -  tETAUt INT) R0SHAP1<1, TNT) *WT  ♦  00*160 

*  ESIG(  IMT)*ODS  U»l(li  1NTJPWTX)  00*170 

EF0RCE13)  -  EF7RCEO)  -  (ETAUtlNT) ♦0SHAP0(2> TNT)  *WT  ♦  CC*180 

*  E$IGUNT)*DDS4AP0t2,INT)RWTl)  00*190 

EFOP.CE  (A)  -  EFORCE  t  A)  -  tETAUt  TNT)  *OSH  API  {Zt  INT)  *WT  ♦  C0A2C0 

♦  ESIGtINT)*DDSHAi>lt2,  INT)*WT1)  00*210 

200  CONTINUE  00*220 

C  00*230 

C  WRITE  ( 6*  2030 )  EFORCEt  1)>  EFOP.CE  ( 2)  t  EFORCE  <  3)/ EFORCE<  *  )  00*2*0 

2030  FQP.MAT 1 5X*  AEI5.6 )  00*250 

C  00*260 

RETURN  00*270 

END  00*280 

SUBROUTINE  S3L''EtA,8,NN,N3*MHAX»KK)  00*290 

DIMENSION  At  1) >  3  1 1 )  .  00*3C0 

00*310 

SOLUTION  OF  SYMMETRIC  BANDED  EQUATIONS  IN  SINGLE  SUBSCRIPT  ARITH.  00*320 

00*330 

MB1-MB-1  00*3*0 

NNN-NM-1  00*350 

IFtKK.EQ.2)  GO  TO  2000  00*360 

IJ-1  00*370 

DO  300  H-l/MNN  00*380 

CC-A(II)  00*390 

IFICC.EQ.O.O)  GO  TO  250  00**00 

41-II+1  00**10 

J2-II+MB1  00**20 

ME-NN-N  00**30 

IF(ME.LT.MBl)  J2-II+NE  00***0 

H«II-1  00**50 

DO  200  4-Jl#J2  00**60 

M-MtMMAX  00**70 

IF(A(J).E0.0.0>  GO  TO  200  00**80 

C«A(J)/CC  00**90 

K-H  00*500 

DO  100  I-J,J2  00*510 

K«K*1  00*520 

100  A(K)«A(K)-C*AtI)  00*530 

A(J)-C  00*5*0 

200  CONTINUE  00*550 

250  CONTINUE  00*5601 

II-II*NMAX  00*570 

300  CONTINUE  00*580 

RETURN  00*590 

2000  11*1  CO*6CO 

DO  500  H-1»UNN  00*610 

CC-A(II)  00*620 

IFtCC.EQ.O.O)  GO  TO  *50  00*630 

4l»II*l  00*6*0 

J2-II*HB1  00*650 

NE-NH-H  00*660 

IFiNE.LT.HBl)J2-II*NE  00*670 

C-B(N)  00*671 
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L»N  004650 

DO  40)  J-J1/J2  0047CO 

L-L+l  004710 

400  B(l>«B(l)-A(J)*C  004720 

3(N)-C/CC  004730 

450  CONTINUE  004740 

II-II+NNAX  004750 

500  CONTINUE  004760 

CC-AIII)  004770 

IF(CC.NE.O.O)  B<NN)«8«NN)/CC  0C4780 

N«HN  004740 

II-H'UX*(HN-2)'t-l  C048C0 

DO  701  I-2/NN  CC4810 

►|«N-1  004820 

IF(A(i:) .EQ.O.O)  GO  TO  630  004830 

J1*II+1  004840 

J2-II+MB1  C04850 

NE-NN-N  004860 

IF  (ME.LT.I1B1)  J2-II+NE  C04870 

C»B(N)  004880 

l»N  004850 

DO  600  J-J1/J2  004000 

L«L*1  004010 

600  C-C-A(4)»B(L)  004020 

B(N)-C  C04030 

650  CONTINUE  004040 

II«II-MMAX  C04050 

700  CONTINUE  004060 

RETURN  C04070 

END  004080 

SUBROUTINE  M0DIFY(N/QjA/8/  NN/  (18/  MM  AX )  004000 

DIMENSION  A(MMAX/NN)/B(MN)  005000 

005010 

MODIFICATION  FOR  PRESCRIBED  ESSENTIAL  BOUNDARY  CONDITIONS.  CC5C20 

005030 

DO  100  J-2/MB  005040 

L»N-J*1  005050 

IF(L.LE.O)  GO  TO  50  005060 

B(L)»B(L)  -  A{ J»  L) *9  C05070 

A<J/L)«0.0  005080 

50  L-N+J-l  005050 

IF(L.GT.UN)  GO  TO  100  005100 

B(L)-B(L1  -  AU/M)*9  005110 

100  AIJ/N)  -  0.0  005120 

B(N)  «  0  005130 

A(l/M)  -  0.0  005140 

RETURN  005150 

END  005160 

SUBROUTINE  GRIO(YC OR D/NODE/ NODES/ NUMEL  )  005170 

DIMENSION  YCGRDI 30 )/ NODE (2/ 30)  005180 

THIS  SUBROUTINE  GENERATES  NOM-OINENSIONAL  CO-ORDINATES  OF  NOOAt  005150 
POINTS.  005200 

Y  -  BAP.  *  Y  /  ( HT/2)  .  005210 

0Y1  «  0.05  005221 

0Y2  ■  0.10375000  005222 

DY  ■  Dv2  ■  005223 

YI  •  -1.0  005230 

DO  100  M  -  1/  NODES  005240 

YCORO(N)  -  YI  005250 

YI  •  YI  ♦  DY  005260 

DY  ■  DY2  005263 

IFCN.GT.-'.AMD.N.LT.IB)  DY  ■  DY1  005264 

100  CONTINUE  C05270 

C  005280 


f 
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:  GENERATE  MODE  NUMBERS 

NO  «  1 

00  200  11  -  1, NUMEL 
N0DE!1,N)  -  NO 
N00£(2>N)  ■  NO  ♦  1 
MO  -  HO  +  1 
200  CONTINUE 
RETURN 
END 

SUBROUTINE  DISTURB <Y CORD, NODE,  EGOOT,  DDDOT,  ALFA, CN» 

♦  £PSILON»NUMEL,NINT,TEMP,tMP,NODES/HTHP  ) 

COMHOH/SHAPGS/ SH APO! 2,  5 )  ,  SHAP1 1 2,  5  ),  OSHAP  .3!  2,  5  )#  OS' 1 A  PI  (2,5), 

♦  00SMAP0(2^51»00SHA|,1(2j5)jSHA(,B(2»5)>0SHAPB(2/5) 

DIMENSION  YCORD!30),NCIDE<2,30),EGDOT(5,30)>DDDOT(5>30) 

♦  ,TMP!31),TEMP!5,30) 

WRITE! 6,2010) 

2010  FORMAT ( 5X, ’MODAL  VALUES  OF  VARIABLES  CAUSED  BY  THE  DISTURBANCE 

♦  5X,  ’ELE.  #  INT.  PT.  #•  ,0X,  •  GAMADOT*, 9X>  •D-D0T’»15X>  ’TEMP.') 
DO  100  MEL  »  1,  NUMEL 

II-  N3DE(1,NEL) 

T2  ■  MODE!  2,  NED 
DO  100  INT  -  1,  MINT 

S  -  SHAPB ( 1* INT) ♦YCORO ! ID  +  SHAPB (2, IMT) *YCOROI 12) 

FAC  -  !1.0-S*S) 

FACN  -  CAC**CM 
EXPV  -  EXP  l-ALFA*S*S ) 

EGDOT!IMT,NEL)  -  0.0 
DDDOT! INT, MEL)  -  0.0 

TEMPIIMTjMEL)  ■  HTHP  ♦  E<>SILOM*FACN*EXPV 

WRITE  16,2000 )NEL , INT, $, EGOOT! INT, NEL), DDDOT l INT, NEL) , TEMP! INT, 
2000  FORMAT! 5X, 215, 4E 15. 5) 

100  CONTINUE 

DO  200  I-l, NODES 
S  -  rCORD!I) 

TMPII)  -  HTMP  +  I ll.-S*S) *+CM) *EPSILON* 1  EXP I-AL  FA*S*S ) ) 

WRITE! 6, 2000)  I, I,S, TMP!I) 

200  CONTINUE 
RETURN 
END 

SUBROUTINE  THERMITEMP,  W0RK,FAC1,  NUMEL, NINT, NTT, DT, TMP 
•*•,  YCORD,  HIDES,  NODS ) 

COMMON/SHAPES/SH AP0!2, 5 ),SHAP1!2,5),DSHAP0!2,5),DSHAP1!2, 5), 
DDSHAP0!2,5),0DSHAP1!2,5)  , SHAPB !  2,  5 ),  DSHAPB!  2,  5  ) 

DIMENSION  TEMP! 5, 30), WORK! 5, 30), NODE! 2,30), YT ( 2), YCORD 130) 
DIMENSION  TMP!31),niSSIP!3),HEAT!2,31), 

♦  THCOMD ! 31, 31), FORCE ! 31 ) 

'  DIMENSION  LD!3)»TH(6),EMAS S (6), WE IGHT ( A ) 

DATA  WEIGHT/O. 347654345137454,  0.  652145154862546, 

♦  0.65 214 5154362546,0 .3 47 854345 1374 54/ 

E  GENERATE  THE  MATRICES 

: 

MB  -  2 

DO  250  I  -  1  ,  NODES 
FORCE! I )  -  0.0 
IFINTT.NE.l)  GO  TO  250 
DO  230  J  -  1,  MB 

230  HEAT! -),T )  -  0.0 

DO  231  J  -  1, MOOES 

231  THCOMD!  J,  I  )  -  0.0 
230  CONTINUE 

DO  ROD  NEL  -  1, NUMEL 
00  30D  I  -  1,  2 


005290 
005300 
005310 
005320 
005330 
005340 
005350 
005340 
005370 
005380 
005390 
0054  CO 
C05410 
0054  20 
CO  5  4  40 
005450 
005460 
•/  C05470 

005480 
005490 
0055C0 
005510 
005520 
CC5520 
005540 
005550 
C05560 
005570 

cosseo 

005590 
NEL) 0056C0 
005610 
005620 
005660 
005661 
005670 
005680 
005710 
005720 
005730 
005740 
0C5741 
CO  57  50 
0C5760 
0057e0 
005790 
005800 
005810 
005820 
005830 
006120 
006130 
006140 
006150 
006160 
006170 
006180 
0C615O 
0062C0 
0062C5 
006210 
006220 
006230 
006  2  40 


» 


II  -  IOOE  ( I,  NEl ) 

vt(i)  -  rcoRoat) 

300  CONTINUE 

IF(NTT .ME. I)  SO  TO  310 
NC  ■  0 

00  400  I  ■  1#  2 
00  400  J  -  1,  I 
HC  -  NC  ♦  1 
THINCJ  •  0.0 
400  ENASStHC)  -  0.0 
NC  -  0 

DO  509  I  -  1  ,  2 
00  500  J  -  1  ,  I 
NC  •  IC  ♦  1 

DJAC  -  (YT!2)  -  YT(1))*0.53 
IFIDJAC.LT. 1.0E-10)  GO  TO  799 
DO  430  IHT  ■  1, HINT 
VT  .  WEIGHT! IMT)*DJAC 
WT1  «  WEIGHT (INT  )/OJAC 

ENASStHC)  •  EHASS(NC)  ♦  SH ARB! I# INT> *SHAP B( J. INT )*WT 
TH(NC)  -  TH(NC)  ♦  OSHAPB(I»IMT)*nSHAPB(J*INT)*WTl 
450  CONTINUE 
500  CONTINUE 
510  CONTINUE 

DO  550  I  •  1*  2 
DlSSIPd)  -  3.0 
OJ AC  -  ( YT (2 )  -  VT (II ) *0. 3  3 
00  550  INT  -  1,  NINT 

0  J  AC  -  DSHAPBd,  INT)*YT(1)  ♦  DSHAP8 ( 2, INT )* YT ( 2) 

VT  -  WEIGHT! INT)  ♦  OJAC 

DlSSIPd)  -  DISSIP!  I )  ♦  SHA<»Bd,INT)*WORKdNT,NEL>*WT 
550  CONTINUE 

C  ASSEM8LE  THE  GLOBAL  MATRICES 
C 

00  600  11-1,2 
I-NODE ( I 1/ NEL ) 

600  LD(II)  -  I 
NC  -  0 

DO  700  II  -  1,2 
I-LD(II) 

FORCE!  I)  -  FORCE  ( I )  ♦  OlSSIPdl) 

IF(NTT.ME.l)  CO  TO  700 
DO  690  JJ  -  1,  II 
NC  »  HC  +  1 
h  -  LO(JJ) 

IR  -  NINOd,  N) 

IC  -  IABSd-D  ♦  1 

HEAT!IC,:R>  -  HE  AT! IC, IR )  ♦  ENASS(NC) 

THCOND(I,H)  •  TUCONO (I,N )  ♦  TH(MC)*FAC1 
THCOND (H, I )  -  THCOND ! I,N ) 

690  CONTINUE 
TOO  CONTINUE 
900  CONTINUE 

00  950  I  ■  1, NODES 
IM1-I-1 

IFd.EO.l)  INI  -  1 
I®1  «  1+1 

IFd.EQ. MOOES)  IP1  •  I 
00  945  J  -  INI, IP1 

94?  FORCEd)  •  FORCEdl  -  THCOND (I,J)*T.HP(J) 

950  COtlTINUE 

IF(NTT.EO.l) 

♦CALL  SOLVE (HEAT, FORCE, NODES, MB, N5, l ) 

CALL  SOLVE (HEAT, =ORCE, NODES, MB, MB, 2) 


006250 
006260 
006270 
006  2  89 
006250 
006  3  CO 
006310 
006320 
006330 
006340 
006350 
006360 
0  C  6  3  70 
006360 
CC64C0 
006410 
006413 
006420 
C06430 
006  4  40 
006430 
C06460 
006470 
006480 
006490 
006 3  CO 
0063C5 
C06510 
006320 
006330 
006340 
006359 
006560 
006390 
0066  CO 
006610 
006620 
006630 
006640 
006630 
006660 
006670 
006680 
006690 
006700 
006710 
006720 
006730 
006740 
006730 
006760 
006770 
006780 
006790 
006792 
006793 
006794 
006793 
006800 
006810 
006812 
006840 
006870 
006880 
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cn  901  -|-1, HIDES  006420 

THP|M)  •  ♦  r  >  .  r  (  M)  *3  T  006930 

900  CONTI  IUE  006960 

C  006  9  SO 

00  o<05  Mr‘.  •  1, -IN-EL  C069 to 

11  •  NODE ! 1,  NE  L I  006970 

12  •  M3DC  I  2.  ‘IEl  )  006494 

00  905  INT  •  1,,|INT  COTOCO 

995  TEHP(INT,NEL >  •  TMP( I  1 )•$  tAPMIl, INT)  ♦  TNP< 12  I *S HAPS < 2 , INT )  0C701O 

•.stjti  0070)0 

4<n  4RITE(6,’4)0)  NEL, MINT, DJAC  007064 

264.)  e)piat(  5*,  ">  n:  *:»  ele.  •>  n, •  iht.  pt.»,  d,  •-•#!  d.)//  cot3;o 

♦  5x,**c  exe of  T.|f  ppigram  is  oeimc  stopped  in  suo.  thorn n co7oeo 

STOP  007070 

ENO  007090 

SUBR  r  *Tr  Mr  03(  V)4)T,  •|I’r»MUHEl»N3DE,  VELD3T,  TGOQT,  ''CORO.  ECO  CT  007  090 

♦  ,  DT,  N2EPQ.  NTT)  007100 

COMMON /SHAPE  I/SH AP0I2, 51  #S<A“U2,5I,0SH6»0(  2,  5 ) .  OSH AR1  <  2,  5  ),  007110 

♦  OOSUAPOCmnoSHAPl^Ol.S  UP8<2,5>,DSHARi(2,S>  007120 

DIMENSION  nooaT< 6. 30), MODE <2, 10 > , V ElOOT ( 11 > , T COOT t 31 ),  007160 

♦  YC39r'(30>,ECI'3T(5»30)  007150 

c  007170 

c  this  suspnriME  cimr-tes  tital  o-iot  and  total  g-oot.  oo7i§o 

c  007190 

C  I F IMTT  ,  L  F . 5 )  t-),2010)  C072CO 

2010  F0RMAT(5<,  'FLc.«  INT.  PT.  »  T TT At  G-0 IT* »5X#* TOTAL  D-DOT»//>  007210 

C  007  2  20 

00  200  NFL  •  1,  Nil 4 EL  C07230 

11  «  NODE ( 1, NEL  )  007231 

12  -  NODE  (  2,  MEL  )  CC7232 

OJACI  •  2.0/  f  vr  j»d  ( 1 2  )  -  TCOROCIin  0C7233 

DO  200  I'|T  •  1,  NIST  007260 

DOOOTI IN’, NEL )  -  004HT ( INT.NEL I  ♦  1 0T*4J AC  1*0 JAC I )•  C073I1 

♦  (ODSHAPm 1, INT) -VELOOTI II  I  ♦  D0SMAP0I2,  INT  )•  VEl  DOT  (  12)  ♦  007)52 

♦  DOS HARK  1, IMT) • TG4HT (II)  ♦  0DSHAP1 ( 2, I  NT )*TCODT ( 1 2 )  I  007323 

EGOOT< INT, NEL)  .  EGOHT ( INT, NEL I  ♦  OTPOJACI*  00T354 

♦  t  0SNAPC(1,INT)«VEL,'0T(T:)  ♦  OSHA»0(2# INT)PVELOOt! 12)  ♦  00735) 

♦  DSHA»1I1,INT) •TGO)T(Il)  ♦  DSMAP1(2,INT)pT6DOT(I2I)  0073)6 

C  IFINTT.LE.))  007360 

C  ♦yRITE(6,3001INEL,I‘i’,EG03T(IflT,MEU,M0D0T(IMT,N€L>  007370 

2000  FORMATOT, 215, 2215.5)  007390 

200  C ONTIMUE  007390 

RETURN  00  76  CO 

ENO  007610 

SUBROUTINE  PLASTIC (FGOaT.DOOOT.PGO'lT.PDOIT,  CNI, GDOT0, $ 16N A,  007620 

♦  CN,A,a,CMU,CN'),0’,TAU0,’EMP,SI0,  NTT, HU NEL, HINT, NORM, TAU  0076)4 

♦  ,  HSI,HCP,HGPOTT,»'-,AMA,EI,PSTAPT#PO)  007631 

DIMENSION  EGDfTM  5, 341, D003TI 5, 30 ) , PGOOT (5 , 3 J> , PDOOT (  ),  30)  007660 

♦  ,  RG  A  U(  5,  100),  SIO,  100),  TEMPI),  )0),TAU(),30),SI  GNA  (  5,  101  CC76  90 

♦  ,W«3RK  <5,301.  »0(  )»  100)  007660 

C  0074  70 

imp$tart.eq.d.o)  go  to  no  00747) 

IFINTT.ME.l)  GO  T0  100  007490 

C  INITIALIZE  S I, PGnGT,POO  )T, TEMP , P GAMA, PO, S IGHA, TAU  AT  EACH  POINT  CC7490 

C  C075CO 

00  SO  NEL  •  1 , MU  ME  L  007510 

DO  SO  HIT  •  l.fllNT  007)20 

SI I IMT , MEL )*HS I  0075)3 

RGOOTI  I'IT,MEL  >  •UGPDT)T  007)60 

PODOT  IINt,NEL)«0.'1  0075)0 

RCANAI  I'IT,NEL).  |C,0  C07570 

RO*riT,MEL  1-0.  3  007)90 

TAU(I*IT,hEL)  •  TAUO  007)90 

5 IGMAl TNT, NFL )  -  3.0  C«f*rO 
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t 


so  comt:  tuc 


0C7610 

C  0C’620 

100  CONTI JUE  007630 

C  0076*0 

DO  503  Hfl  •  l.MUHEL  007630 

00  300  IMT  .  l,MtMT  007660 

SII-StlUtT.HELI  007670 

TNP.TEWPIINT.NEU  0076M 

CKA*A1  -  Cll.O*Stt/SIO)**CH»  007630 

CK  APA  •  CKAPA1  •  (X.O  -  A*THP)  0077CO 

ESTHEl  •  SOOT  (  TAU(  IMT,  NEU*TAW«  IHT.NEL  >  ♦  CC7710 

*  S:G*Ut:HT,MEL>»SIGNACIMT,»lEL)  )  007711 

OT  ■  ESTPES/CKAPA  0C7720 

tFM’.L’.O.OI  07  73  4,6!  007723 

:F(a7.LT.1.0)  G7  7n  *33  0C7730 

*2  •  OT**CHI  -1.0  007730 

GAIA  •  H2/ t3PESTHES*G03T0)  CC7760 

S1O07  •  GAflA  •  EST6ES*ES7*ES/CKAPA1  007770 

GO  TO  *60  CC7780 

630  CONTINUE  007763 

C  WHITE ( 6,2100)  IMT, NEL, 3T  007760 

2100  FOPHATJJX,  '7M£  INT.  PT.<,!5,'  OF  ELE.  A  '.n.'JS  UNLOADING.1/  C073C0 

♦  5X,»T4I>/P.APA  •  •*  El 5 ,5)  CCT810 

gana  •  o.o  007*20 

SIOOT  •  7.0  007830 

*60  CONTIN'IC  007360 

GP007  •  GAIA  •  TAUCNT.MEL)  007830 

7»0JT  .  GANA  •  S ’G*i( INT*NEL 1  007860 

PGOOTI IN’.NEL )-G»OOT  007870 

•HOOT!  IN’,  NEL  I  -  08077  007880 

SKI*)’, MED  •  SIIINT.HED  ♦  SIOQTpO’  007890 

PGANA! IN’, NED  •  PGANAI IM’.NEL )  ♦  G»OOT»D’  007900 

PO(INT.NEL)  •  PO ( TNT*  NEL I  ♦  0»03T*0T  CC7910 

C  WH I TE  (  6, 2000 )  NEL* IN’. SIOOT, GPOOT,OPOOT  007920 

2000  F7«NAT(3X*2I3,3E13.3)  007930 

WOHKIINT.NEL)  •  CKAPA1  •  SIOOT  007960 

TAUOOT  .  CNUPIEC70T(  INT.NEU  -  GPOTT )  007930 

SIGOOT  •  CMU* (DOOQTI TNT, NEL )  -  0»D0T)  007960 

TAU(  INT,  MEL )  •  TAIKIHT.NEL)  ♦  TAUD0TP7T  007970 

3 IGNA( IN’, MEL )  •  SIGNAI TNT.NEL )  ♦  SIOOOTPOT  007980 

300  CONTINUE  007990 

GO  77  *62  007991 

661  WHITE ( 6, 2010 )  OT  007992 

2010  FQHNA7 ( /// /3X, 1 TNE  EXECUTION  OF  THE  PHOGHAN  IS  BEING  007997 

•  S’JPPEO.  QT  •  ' , E 13.3 )  00’99* 

S’OP  007995 

*62  CONTINUE  007996 

HE’UHN  00 80 CO 

EMO  008010 
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